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Abstract. 

Conventional tomographic techniques are becoming increasingly infeasible for recon- 
structing the operators of quantum devices of growing sophistication. We describe a novel 
tomographic procedure using coherent states which begins by reconstructing the diagonals 
of the operator, and then each successive off-diagonal in a recursive manner. Each recur- 
sion is considerably more efficient than reconstructing the operator in its entirety, and each 
successive recursion involves fewer parameters. We apply our technique to reconstruct the 
positive-operator-valued measure (POVM) corresponding to a recently developed coherent 
optical detector with phase sensitivity and number resolution. We discuss the effect of various 
parameters on the reconstruction accuracy. The results show the efficiency of the method and 
its robustness to experimental noise. 
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1. Introduction 

Quantum detectors inform our classical world of the underlying quantum world through a 
set of operators known as positive-operator-valued measure (POVM). In practice, the success 
of many quantum applications rely on certain knowledge of measurement apparatuses [1, 

2, 3]. Successful applications of sophisticated detectors rely on a complete and accurate 
knowledge of the detector, i.e. detector characterization. Detector characterization can 
be implemented in two different ways. One is synthetic, wherein each constituent of a 
detector is carefully calibrated before being incorporated into a sophisticated physical model 
of the measurement process. As quantum technologies evolve into increasingly complicated 
systems, so do quantum detectors, which makes synthetic characterization progressively less 
feasible. Additionally, any coupling with external degrees of freedom not incorporated into 
the theoretical model may make the characterization fail [4, 5]. A fundamentally different 
approach is taken in quantum detector tomography (QDT) [6, 7, 8, 9], where the unknown 
specifics of a detector are characterized in a largely assumption-free way: here, the POVM of 
a detector are reconstructed from the outcome statistics in response to a set of tomographically 
complete certified input states. 

To date QDT has been successfully applied to avalanche photodiode (APD) [ ], time- 
multiplexed photon-number-resolving detector (TMD) [9, 10, 1 1], transition edge sensors [12] 
and superconducting nanowire detectors [13]. These detectors are phase-insensitive, i.e. 
they can only measure the mixture of the photon-number states, not the coherence among 
them. Accordingly, the POVM of these detectors have only non-zero matrix elements on 
the main diagonals, and the number of parameters to be decided is proportional to d. Here 
d is the dimensionality of the Hilbert space, and for optical detectors can be estimated as 
the number of photons to saturate the detector. Yet a large number of effects characteristic 
for quantum mechanics, including entanglement, violation of local realism [ ], measuring 
non-classical correlations of radiation fields [15], test of macroscopic realism [16] etc., relies 
on quantum coherence. The effort to harness and exploit quantum coherence brings the 
prosperity of quantum information processing and quantum metrology. Moreover, exploration 
and utilization of the full Hilbert space of a quantum system requires a detector capable 
of implementing a tomographically complete set of measurements [ ]. Therefore optical 
detectors that can access quantum coherence among photon-number states, i.e. phase- 
sensitive detectors, for example strong- and weak- field homodyne detectors [ ], are crucial 
not only for quantum applications, but also for test of fundamental theories of quantum 
mechanics. Phase- sensitivity comes with the non-zero off-diagonal matrix elements of the 
POVM. Thus tomography of a phase- sensitive detector requires the estimation of number 
of parameters proportional to d 2 . For practical detectors, d can range from 10 2 to 10 5 , and 
d 2 from 10 4 to 10 10 . For example, a weak-field homodyne TMD with 9 time bins requires 
1.8 x 10 6 parameters to completely describe its POVM [5], which is about two orders of 
higher than the largest tomography that had been performed until then [19]. Such large set 
of parameters represents a considerable challenge to the characterization of phase-sensitive 
detectors. 
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In this work we explore potential solutions to QDT of phase-sensitive quantum detectors. 
In particular we introduce an algorithm that allows to reconstruct the POVM recursively, with 
no more than d parameters per recursion. Simulations with the QDT of weak-field detectors 
demonstrate the robustness of this algorithm. 

2. Definition of the problem 

QDT is performed by preparing a set of known probe states {p m } incident on a quantum 
detector and observing the detector outcomes. The probability of registering outcome n is 
given by the Born rule 

Vn\m = tr(p m Il n ), (1) 

where {II n } is the POVM of the detector with n = 0, . . . , N — 1, and N is the number of 
possible outcomes of the detector. In practice the experiment is repeated for each of many 
identical copies of the probe states, and the frequency f n \ m for each measurement outcome 
n occurring when probe state p m is used is recorded. Then p n \ m can be estimated from the 
relative frequency p n \ m = f n \m/ J2 n fn\m- One can then invert Eq. (1) to find Il n . For a finite 
number of repetitions, there are always fluctuations in the estimation of p n \ m , therefore the 
inversion should normally be preformed with convex optimization. 

A key requirement is that the set of probe states must be tomo graphically complete. 
However, it is also important that the set of probe states are experimentally feasible. That 
means that the states should themselves be well-characterized, and that a large variety should 
be available with high precision. There are proposed methods to generate the probe states 
through quantum correlations [8, 20]. Yet with current quantum optical sources it is very 
hard to generate the probe states strong enough to saturate the detector under test. For 
photodetector measurements, there is a more straightforward option. The set of coherent 
state vectors \a) of an optical beam are ideal candidates, where a = \a\e td is the complex 
amplitude. They are overcomplete in the sense that two different coherent states are not 
orthogonal with each other yet any quantum state can be decomposed on the set of coherent 
states. Therefore coherent states can form a tomographically complete set by transforming 
their amplitude (by means of optical attenuation) and their phase (with a simple delay line). 
Importantly, they are generated very easily by a laser. 

With coherent states as input, the probabilities are given by 

p n \ a = (a;|n n |a) = 7iQ n (a), (2) 

where Q n (.) is the Q-function of the detector POVM elements n„. This is equal to the 
Husimi representation of the POVM, and is uniquely and invertibly related to the POVM. 
To reconstruct U n one can write both |a) and IT n in the photon-number basis and truncate the 
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expansion at d — 1, where d — 1 is the number of photons that saturate the detector. 

|«) = e"H 2 / 2 X:^b), (3) 

3=0 V J 

d-1 

n«= J2< k \j)(k\- (4) 

j,k=0 

Then Eq. (2) can be written as 

d—l I \j-\-k 

= e " |Q ' 2 E %r^ ik ~ ])d ^- (5) 



/JIM 

We can relabel Eq. (5) in s = kd + j ' + 1 (1 < s < d 2 ), with j = [(s — 1) mod d] and 
^ = ( s ~ J ' ~ 1)/^- F° r M probe states, there are M x N linear equations, which can be 
written in a matrix form 

P = FTl (6) 

where P is an M x iV matrix with elements p m > n = p n \ am , F is an M x <i 2 matrix with 
elements 

In/ |i( s )+^( s ) 

and II is a <i 2 x iV matrix with elements 7r s,n = 7iv^' fe . In practice where the experimental 
noise is taken into account, the POVM set can be estimated from Eq. (6) with convex 
optimization subject to the constraints 

n n > o, (8) 

N-l 

J2 U n = I, (9) 

n=0 

where / is the identity operator. One common approach is the least square estimation 

min||P-Fn|| 2 , (10) 



where \\A\\ 2 = ^Tr(A^A) is the Frobenius norm. The reconstruction problem effectively 
deconvolves a coherent state from the statistics to obtain the POVM set. This is an ill- 
conditioned problem, as seen by the large ratio between the largest and smallest singular 
values of the matrix F. This makes the POVM extremely vulnerable to small fluctuations 
in the measurement statistics. Such instability can be overcome by adding a regularization 
function g(TL) to the optimization [9, 10], therefore the problem is modified as 

min{||P-Fn|| 2 + 2(n)}, 



N-l 



subject to il n > 0, ^2 n ™ = L C 11 ) 



n=0 
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For a phase-insensitive detector with finite detection efficiency, one would expect the variation 
of the diagonal matrix elements to be smooth, therefore a regularization function known as 
Tikhonov regularization [ ] is applied 



This limits the variation between adjacent elements along the diagonal matrix elements. Yet 
for a phase- sensitive detector a regular function is not easy to find: even as each of the leading 
diagonals are smooth, the relation among different leading diagonals can be arbitrary. 

An alternative approach for convex optimization is maximum likelihood estimation, 
which was also proposed for QDT [7]. Maximum likelihood alleviates the requirement of the 
regularization function. However, its convergence speed is normally not high. Moreover, both 
the maximum likelihood estimation and the least square estimation in Eq. (11) requires the 
reconstruction of the whole POVM matrices at the same time. When the size of the matrices 
becomes large, the problem becomes infeasible. For example, the estimation of a POVM set 
with 9 elements each of which is a 50 by 50 matrix is already a hard problem for the capability 
of current multi-processor desktops (2xQuad Core 3GHz, 8GB RAM). 

The engineering of large entangled quantum states and development of sophisticated 
quantum operations has set a challenge for standard quantum tomography techniques. There 
has been increased interest in the development of novel algorithm with improved efficiency 
for special situations. In particular, there are process tomography schemes that allows to 
selectively reconstruct the state or process matrix partially in each run. Several of them use a- 
priori knowledge about the state such as their symmetry [22, 23, 24], or simply reconstructing 
a subset of the entire state or process [25, 26]. Using improved techniques from classical 
signal processing have also become common, such as compressed sensing [27, 28]. In the 
following, we introduce a novel algorithm that reconstructs the POVM elements recursively 
in multiple runs. 

3. The detector model 

The algorithm discussed in this work can be universally applied to the tomography of any 
quantum detector. To illustrate its working in this work we consider a simple example of 
phase- sensitive detector, the weak-field homodyne APD. The configuration of such detectors 
together with a schematic tomography setup is given in Fig. (1). The probe state is prepared 
by the phase modulation and amplitude modulation of the output of a laser system. The weak- 
field homodyne detector is shown in the black box, where the input state interferes with a local 
oscillator (LO) and detected by a photon-number-resolving or photon-counting detector. For 
a weak-field homodyne APD, there are two possible measurement outcomes, no-click and 




(12) 
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click events, and the corresponding POVM elements il and IIi are given as [18]. 

e -KI 2 



c=0 d=0 

x (a* + a f ) c K - a t ) d |0)(0|(a L + a) c (a L - a) d , 

iii = / — n . 



(13) 
(14) 



where a L is the complex amplitude of the LO and t/apd is the detection efficiency of APD. 
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Figure 1. The configuration of a weak-field homodyne detector and its tomography setup. 
A set of probe states are prepared by the phase-modulation (PM) and amplitude-modulation 
(AM) of the output of a laser. The magnitude of the probe state is adjusted by a half-wave 
plate (HWP) followed by a polarizing beam-splitter (PBS) and neutral density filters (ND). 
The phase of the probe state is controlled by a piezo translator. The setup of the weak-field 
homodyne detector (WHD) is shown in the black box. 



4. A selective algorithm with Glauber-Sudarshan P-function 

Before we proceed to the recursive algorithm, we consider a more straightforward selective 
algorithm. Each matrix element 7r^ fc of U n is given by 

<< fc =Tr(|A;>(j|n n ). (15) 

Using the Glauber-Sudarshan decomposition of \k)(j\ 

\k)(j\=2 J P^ k (a)\a)(a\d 2 a, (16) 

we have 

< - 2 /^(„)( t , W « ) A = 2 I / P »«, H A. ,17) 

In principle we can estimate each individual matrix element separately with either 
the exact form of Pi' k (.) which contains the derivative of Dirac-delta function or the 
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approximated regular form [29]. Similar method has been used for quantum process 
tomography [30, 31, 32]. However, due to the singularity of P-function, this scheme is 
extremely sensitive to the noise in the measured Q-function of the POVM element, rendering 
it infeasible for practical QDT. 

As an example, we consider the reconstruction of no-click POVM of a weak-homodyne 
APD with the reflectivity of the beam-splitter of 0.5, LO intensity |ctLo| 2 = 5 and quantum 
efficiency of the APD 60% (overall detection efficiency 30%). We choose the probe states 
that sample the phase space from X, Y = —10 to X, Y = 10 with a step size 0.05, where 
X and Y are the two quadratures of an optical field. We assume that for each probe state 
we run the experiment / times, then the expected frequency to get the no-click event is 
(fo\a) — n Qo(oi)f. In practice there are many experimental imperfections that may induce 
fluctuations to the measurement results. In this work we only consider the most fundamental 
fluctuation due to the random nature of the outcome of each measurement process, and 
simulate it by assuming that f \ a is a random number with a binomial distribution, and 
assigning the experimentally measured Q-function as Q exp (a) = fo\ a / f- 

The results are shown in Fig. (2). Without experimental fluctuations, the reconstructed 
POVM matches almost perfectly with the theoretical prediction. Yet when there is the 
presence of experimental noise, the results deviates from the theoretical prediction very 
quickly: for / = 10 5 as in Fig. (2(c)), the reconstructed POVM element is not even a physical 
measurement operator. Only when the number of measurement is large enough / = 10 10 
and thus the experimental fluctuations is small, the reconstructed POVM element is close to 
the real one. The results presented here is reconstructed up to the 4-photon component. The 
P-functions of higher photon-components are more singular, and the reconstruction is more 
sensitive to experimental fluctuations. Therefore, this method has a serious problem for its 
scalability. For practical QDT, we need to seek another solution. 

5. Recursive reconstruction of the POVM set 

5.1. Outline of the recursive reconstruction 

In this section, we discuss a novel recursive method for the tomographic reconstruction of 
quantum operators. We begin with Eq. (5). Before relabeling, we first integrate over the 
probe state phase 9. With 



The left side of Eq. (19) is a partial integration of the experimental results, while the right side 
involves only the main diagonals of the POVM. Eq. (19) can be interpreted as using phase- 
randomized coherent states as input to the detector. Since the input states are completely 
mixed, the measurement process only involve the main diagonals of the POVM. In a practical 




(18) 



we have 




d-l 



(19) 
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experiment, one should change the integration on the left side to summation. The probe 
states should be prepared with M a different amplitudes and for each amplitude there are M p 
different phases. In total there are M = M a M p probe states with the complex amplitudes 
oc u ,v — |a«|e l6,u ", with u = 1, . . . , N a and v = 1, . . . , N p . Therefore the integration on the left 
side of Eq. (19) can be approximated as 




For reconstructing the off-diagonals we first note that POVM elements are Hermitian and it 
is sufficient to reconstruct just the upper or lower off-diagonals. We multiply Eq. (5) by e~ lW 
and integrate over 9. Since 

/ \^ e d9 = 2nS kJ , (21) 
Jo 
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we have 

„- / Pn\«e- iie dO = r e >' a . |Qt| 7T^. (22) 
Again for practical experiment the integration should be substituted with a summation 



- / Pn\ a e- iie d0 « — ^^.e-^-. (23) 



with an error 



27T 3 rf 2 (pn| a e-™) 27T 3 / 2 

3M 2 ~ 3M 2 ' 1 j 

Eq. (22) includes the situation in Eq. (19) when I = 0. For each /, there are M a equations. 
As has been done in Eq. (6), we can write them in a matrix form pW = p( l )f[W^ with 
pW an M a x N matrix, PW an M a x d matrix, a d x N matrix, and the coefficients 
given by Eq. (22). Comparing with Eq. (6), all the matrices involved here are significantly 
smaller. With the presence of the experimental fluctuations, the reconstruction becomes a 
convex optimization problem, in fact a semi-definite problem, 

min{||P w - P (i) n w || 2 + g(fl {l) )}, 

N-l 

subject to Ii n > 0, n„ = I. (25) 

n=0 

Since this is a convex optimization problem, there is only one minimum value which can be 
found with YALMIP toolbox for Matlab [33] with the solver SeDuMi [ ] utilizing primal- 
dual interior point methods [ ]. For the examples discussed in this paper, the calculation 
converges to its minimum in less than 40 iterations which takes about 1 seconds on a multi- 
process desktop (Dual Core 2GHz, 2GB RAM). This allows us to reconstruct the POVM 
recursively for I = 0, . . . , d. For I = 0, the second condition is that the summation of the main 
diagonals of all the POVM elements equals to 1, while for I ^ 0, this condition means that 
the summation of the Zth leading diagonals of all the POVM elements equals 0. The positivity 
condition should be enforced recursively based on Sylvester's criterion, which states that a 
matrix is positive if and only if all of its principal minors are positive. For I = this requires 
all the matrix elements on the main diagonals to be positive. Now we derive the condition for 
/ > 0. We start with / = 1. In Eq. (26) we show the matrix II n where the diagonal elements 
(green) have been determined using Eq. (19) and the first row of off-diagonals (red) is to be 
determined with Eq. (23). Any other entry is unknown and will not be reconstructed at this 
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(26) 



For an input state vector of the form |$) = a \j) + b\j + 1) the effective submatrix of 
n n is given by 



IE' 1 



7T; 



3,3 



71 



■3,3+1 



Ti 



(27) 



which needs to be positive. We marked these submatrices with blue, orange, and black 
lines for the j = 0, j = 1, and j = M — 1 cases in Eq. (26). Thus we imposed the 
additional constraint that all diagonally centered 2x2 submatrices of IT^ 1 need to be positive 
for the reconstruction of the first off-diagonal. This condition is satisfied if and only if the 
determinant det(n{' 1 ) is positive, which implies 



IT; 



^' +1 | 2 <<'X + W VnJ. 



(28) 



For the following reconstruction of the /th leading diagonal we impose a similar 
constraint on the (1 + I) x (1 + I) submatrices start with ix™ , illustrated in Eq. (29) 



TT' 



■3,3 



7T: 



3,3+1 



\ 



(29) 



Ti 



3,3+1* ... ... n j+l,j+l 



where only tt^^ 1 is unknown. It is required that IT^ is a positive matrix. Since 
constraints in the previous steps ensure all its principle minors are positive, this condition 
is equivalent to that its determinant det(n{ ,z ) is positive. To facilitate numerical calculations 
we derive the bounds on tt^^ 1 , which can be done by noticing that det(n{' / ) is a quadratic 
polynomial of it^ +l 



det(n^) = A x (tt^'+O 2 + B x tt^ +1 + C. 



(30) 
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It is easy to see that A is positive, since A is the product of the elements along the anti-diagonal 
and is Hermitian. Therefore det(II^) > implies that 



2A - n - 2A 

The value of A, £> and C can be easily estimated from Eq. (30) by substituting = ±1,0 
into and calculating the determinant numerically. 

5.2. The number of leading diagonals I 

To reconstruct the full POVM matrices, we should run the calculation in Eq. (25) until 
I — d — 1. As can be seen from Eq. (24), for higher I it requires increased number of phases 
M p to reduce the numerical error. Yet, in practice the number of leading diagonals can be 
estimated during the calculation. From the positivity condition, one has 



TP 



3 ' j+ f <TT 3 ' j 7r J n +l ' j+l . (32) 



Therefore, after the reconstruction of the principle diagonals, we can put a bound on the 
number of leading diagonals to be reconstructed. Moreover, in any practical detector there 
is always a finite fluctuation of the reference phase (with a fluctuation of 2n for a phase- 
insensitive detector), which will further reduce the number of leading diagonals, as shown 
below. In fact, this phase noise will ensure that the entries of the POVM elements decay 
exponentially away from their main diagonal. 

Assume the reference phase has a Gaussian distribution with width of 5 > 0. Instead of 
having a POVM IT n we have 

W n = —5= dZR(O j n n R(0 exp(-£ 2 /(25 2 )), (33) 

where R(£) is the rotation operator in phase space with angle £. The matrix elements of W n 
are given by 

n 



6y2n J- n 



5V2n Jo 

' 1 > i 



<5V2^ 
T 3d+i 



dt(j\R0n n R(£) \j + 1} exp(-e/(25 2 )) 
dm^n\j + l)exp(-e/(25 2 )+UO 
^ d^exp(-e/(25 2 )+UO 

J —n 



= ^= f d£exp(-e/(25 2 )) «*(/£). (34) 
oy 111 J- n 

Intuitively, if the fluctuation of the phase reference is small, i.e., 5 <C tt, the last integration in 
Eq. (34) can be approximated as 

^exp(-£7(2<5 2 ))cos(ZO « / ^exp(-£7(2<5 2 )) cos(Z£) = 5v / 2^exp(-/ 2 5 2 /2). 

-7T J — OO 

(35) 
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The intuition of exponentially decaying coefficients can be made rigorous as follows. One 
has, for w.l.o.g. I even, 



^exp(-£7(25 2 ))cos(ZO 



< 5V27rexp(-ZV/2) 



+ 2 



^exp(-£7(25 2 ))cos(ZO 



5V27rexp(-Z 2 cr/2) 



^exp(-(£ + 7r)7(25 2 ))cos(^) 



< (5V27rexp(-/ 2 (572) 



+ 2 



^exp(-e 2 /(25 2 ))cos(/0 



= 25V2^exp(-/V/2). 

Thus, the matrix elements of satisfy 

\7r'i> j+l \ <2|vr^ +i |exp(-/ 2 5 2 /2). 



(36) 



(37) 



The Zth leading diagonal is decreased by a factor of 2 exp(— 1 2 5 2 /2). With the increase of 
I this factor increases therefore reduces the number of significant leading diagonals in IT^, 
leading to I <C d. That is to say, the effort of reconstruction up to a constant error is of order 
0(d) instead of 0(d 2 ). For example, with a phase fluctuation of 10 degrees the 18th leading 
diagonal is reduced to 1% of that with no LO phase fluctuation. 

Another reason for the reduction of the required calculation for the leading diagonals 
comes from one of the major point of performing detector tomography: to predict the response 
of the detector with various input quantum states. For situations involving input states with a 
fixed photon number N, like NOON states [35] or Holland-Burnett states [ ], we only require 
iV leading diagonals of the POVM elements to predict all measurement outcomes. Due to the 
lack of bright quantum sources, iV is usually small (less than 8). 



5.3. Regularization 

The numerical stability of a reconstruction algorithm is one of its vital certificates. Numerical 
instability has been a common problem in tomography [36, 37], particularly so in using phase 
space data from homodyne tomography to reconstruct operators in the Fock spaces [38]. Tools 
such as pattern functions [39, 40, 41] exist that can bridge this gap. They are however, hard 
to identify and cumbersome to work with [42]. The use of maximum likelihood functions has 
also been suggested for detector tomography [7, 8]. Unfortunately, as mentioned earlier, the 
speed of the convergence of such algorithms is not generally guaranteed to be high, becoming 
exponentially slow for certain problems. 

We strike a balance by developing a recursive algorithm that is efficient by virtue of being 
cast as a semi-definite programme, as is evident from the convex function to be minimized, 
and the linear constraints in Eq. (11). Unfortunately, this still leaves us with an ill-conditioned 
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problem, primarily due to extremely large ratio between the largest and the smallest singular 
values of the matrix F®. This is a consequence of the large range of coherent state amplitudes 
needed to cover the entire dynamical range of the detector in the Fock space. The most 
common outcome of this ill-conditioning is to result in reconstructed POVMs that have sharp 
discontinuities [10]. As shown in Eq. (12) this can be resolved by a smoothing function 
or Tikhonov regularization [21]. We will next discuss how this mathematical technique is 
physically enforced in realistic detectors. 

Most realistic optical detectors have finite efficiencies which enforces a certain degree 
of smoothness in their corresponding POVM representations. If a lossy optical detector has 
a POVM element with non-zero amplitude \m)(n\ it will also have a non-zero amplitude in 
\m + l)(n + 1|, |ra + 2)(n + 2\, . . . ,\m + K)(n + K\, decreasing with K. If the detector has 
a finite efficiency rj, it will impose some smoothness on the distribution 7r^ ,fe . That is because 
if G(k) is the probability of registering k photons and H(k') is the probability that k! were 
present, then the loss process will impose 

G ^ = E (T) ^ - vf- k H(k'). (38) 

This motivates an immediate generalization of Eq. (12) to that in Eq. (25) as 

g(U^)= 1 J2\< J+l -< +hJ+l+1 \ 2 - 09) 

While 7 is a free parameter introduced into the problem for numerical smoothness, we 
show that the outcomes of our reconstruction procedure are fairly insensitive to the actual 
value of the parameter. Fig. (3) presents the effect of the regularization condition for the 
reconstruction of the no-click POVM of a weak-homodyne APD with the reflectivity of the 
beam-splitter of 0.5, LO intensity | ollo 1 2 = 5 and quantum efficiency of the APD 60% (overall 
detection efficieny 30%). We vary the weight of the regularization condition for two orders of 
magnitude. In addition to the fidelity, we also calculate the relative error of the reconstructed 
POVM | |n^ ec - n* he | 1 2 / 1 |n* he | | 2 . The results are presented in Table (1), which show that the 
change of the reconstructed POVM elements due to the change of regularization strength is 
small. This confirms that the main effect of the regularization condition is to suppress the 
ill-conditioning and noise while leaving the POVM fitting unaffected. 



Table 1. Sensitivity of the reconstruction procedure to the choice of parameter 7. No-click 
event of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity 
|qlo| 2 = 5 and quantum efficiency of the APD 60% (overall detection efficieny 30%). 
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As a comparison, we also calculated the reconstruction of the no-click POVM of a weak- 
homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity |«lo| 2 = 5 and 
quantum efficiency of the APD 20% (overall detection efficieny 10%) and that of a weak- 
homodyne APD with the reflectivity of the beam-splitter of 0.1, LO intensity |«lo| 2 = 5 
and quantum efficiency of the APD 90% (overall detection efficieny 81%). The results are 
shown in Figs. (4) and (5). Calculated fidelities and relative errors are presented in Tables (2) 
and (3). We can see that regularization works very well for moderate and low detection 
efficiencies, while its performance decreases if the detection efficiency is very high since the 
corresponding POVM elements are not smooth any more. On the other hand, one can infer 
the detection efficiency from the differences between the reconstructed results with different 
regularization strengths. If such difference is large, one should utilize a reduced regularization 
strength in the reconstruction. 



Table 2. Sensitivity of the reconstruction procedure to the choice of parameter 7. No-click 
event of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.5, LO intensity 
|a_Lo| 2 = 5 and quantum efficiency of the APD 20% (overall detection efficieny 10%). 



7 


Fidelity 


Relative error 


0.1 


99.85% 


1.34% 


1 


99.87% 


1.28% 


10 


98.87% 


2.32% 




Table 3. Sensitivity of the reconstruction procedure to the choice of parameter 7. No-click 




event of a weak-homodyne APD with the reflectivity of the beam-splitter of 0.1, LO intensity 




|a£o| 2 = 5 and quantum efficiency of the APD 90% (overall detection efficieny 81%). 


7 


Fidelity 


Relative error 


0.1 


99.87% 


1.82% 


1 


96.95% 


8.29% 


10 


71.08% 


45.96% 



5.4. Reconstruction of the POVM of a weak-field homodyne APD 

To discuss the performance of the recursive reconstruction method, we numerically simulate 
the reconstruction of the POVM set of a weak-homodyne APD with the reflectivity of the 
beam-splitter of 0.5, LO intensity |ctLo| 2 = 5 and quantum efficiency of the APD 60%. We 
choose the intensity of the probe state \a u \ 2 from to 100 photons with a step size of 0.5 
photon. This is sufficient to saturate the detector response. For each intensity we consider 
probe phases distributed uniformly between and 2n, i.e. 6 U V = {0, 2w/M p , . . . , 2(M P — 
l)n/M p }. Again we only consider the fluctuation induced by the random nature of the 
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measurement process. We assume that for each probe state we run the experiment / times 
and simulate the experimental noise by assume / i Q is a random number with a binomial 
distribution. In Fig. (6) we show the theoretical prediction and the reconstructed POVM for 
the no-click event with M p = 5, 20, 40 and / = 10 5 . To illustrate the results we show each 
leading diagonal separately up to I = 3. The reconstruction is done up to 150 photon-number 
component and is only displayed to 25 photon-number component for clarity. We calculate 
the fidelity between the reconstructed POVM I1q cc and theoretical prediction U)^ e 



which are 87.04%, 98.19% and 98.32% for M p = 5, 20 and 40 respectively. The change in 
fidelity can be further elucidated by the red bars on top of the reconstructed POVM element 
which indicate the distance from the theoretical prediction. From the results we can see that 
although all three phase settings give almost the same results for the principle diagonal, for 
higher I it requires more probe phases for an accurate reconstruction. This is due to the 
numerical error for the calculation of the integral given in Eq. (24). This on the other hand 
shows a practical advantage of the recursive QDT The probe phase setting can be decided by 
the elements in the POVM matrices to be reconstructed. If we are only interested in the low 
leading diagonals, we can greatly reduce the number of probe phases from that needed for a 
complete reconstruction of the POVM. 

In Fig. (7) we show the performance of the recursive QDT under different level of 
experimental fluctuations. Here for each probe intensity we adjust M p = 40 phases. The 
method discussed in Sec. 4 requires / = 10 10 to achieve a satisfactory accuracy. As a 
comparison, the recursive QDT is very robust against the experimental fluctuations: a decent 
accuracy can already be achieved for / = 10 3 (fidelity with the theoretical prediction 98.27%), 
with further improvement for / = 10 5 (fidelity 98.32%). Depending on the repetition rate of 
the detector and the laser system for LO and probe state, this requires only several millisecond 
to one second for each probe state. 

6. Conclusion 

Phase-sensitive quantum-optical detectors are crucial to fully exploit the fundamental features 
of quantum physics and to optimally utilize optical telecommunications channels [43, 44, 45]. 
The success of these applications relies on the accurate knowledge of detectors. Yet as 
quantum-optical detectors become more sophisticated, normal parameters like detectivity, 
spectral sensitivity and noise-equivalent power are not sufficient to provide a complete 
specification of the detector. Moreover the complex structures of detectors and the coupling 
with external degrees of freedom make the conventional characterization of these detectors 
less feasible. Quantum detector tomography, a black-box or device-independent approach 
for the complete characterization of quantum detectors, provides a universal solution to 
this problem. Full characterization enables more flexible design and use of detectors, be 
they noisy, nonlinear, inefficient or operating outside their normal range. However, the 




(40) 
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large number of parameters associated with the tomography of coherent quantum detectors 
presents a technical challenge. This challenge is becoming increasingly typical as quantum 
devices grow in sophistication. In this work we present a novel recursive reconstruction 
algorithm to overcome this problem. Aided by numerical simulations, we have demonstrated 
successful reconstructions of the POVM of a weak-field homodyne APD. The results show 
the flexibility of the algorithm and its robustness to experimental noise. The capability to fully 
characterize coherent quantum-optical detectors paves the way to study genuine quantum 
features, including wave-particle duality, super- sensitivity etc., of a measurement process. 
It allows the benchmarking of the performance of quantum-optical detectors for various 
quantum applications and sheds new light on the assessment and verification of more complex 
detectors. We also hope that recursive quantum tomography provides an efficient procedure 
for quantum tomography in other quantum state and process characterization problems. 
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Figure 3. Reconstructed POVM element for the no-click event of a weak-field homodyne 
APD (with the reflectivity of the beam-splitter of 0.5, LO intensity |a_Lo| 2 = 5 and quantum 
efficiency of the APD 60%) under different level of regularization 7 = 0.1, 1, 10. The 
simulation is done for M p = 40 and / = 10 5 . We demonstrate each leading diagonal 
separately up to / = 3. Red bars on top of the reconstructed POVM element indicate the 
distance from the theoretical prediction. 
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Figure 4. Reconstructed POVM element for the no-click event of a weak-field homodyne 
APD (with the reflectivity of the beam-splitter of 0.5, LO intensity |a_Lo| 2 = 5 and quantum 
efficiency of the APD 20%) under different level of regularization 7 = 0.1,1,10. The 
simulation is done for M p = 40 and / = 10 5 . We demonstrate each leading diagonal 
separately up to / = 3. Red bars on top of the reconstructed POVM element indicate the 
distance from the theoretical prediction. 
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Figure 5. Reconstructed POVM element for the no-click event of a weak-field homodyne 
APD (with the reflectivity of the beam-splitter of 0.1, LO intensity |a_Lo| 2 = 5 and quantum 
efficiency of the APD 90%) under different level of regularization 7 = 0.1,1,10. The 
simulation is done for M p — 40 and / = 10 5 . We demonstrate each leading diagonal 
separately up to I = 3. Red bars on top of the reconstructed POVM element indicate the 
distance from the theoretical prediction. 
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Figure 6. Theoretical prediction and reconstructed POVM element for the no-click event of a 
weak-field homodyne APD. We consider three different probe phase settings M p = 5, 20, 40. 
For each probe state we assume the experiment is run / = 10 5 times, and simulate the 
experimental fluctuation with a binomial distribution. We demonstrate each leading diagonal 
separately up to the I = 3. Red bars on top of the reconstructed POVM element indicate the 
distance from the theoretical prediction. The results are constructed with the weight of the 
regularization function 7=1. 
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Figure 7. Reconstructed POVM element for the no-click event of a weak-field homodyne 
APD under different level of experimental fluctuations. For each probe state we assume the 
experiment is run / = 10 3 and 10 5 times, and simulate the experimental fluctuation with 
a binomial distribution. We demonstrate each leading diagonal separately up to the I = 3. 
Red bars on top of the reconstructed POVM element indicate the distance from the theoretical 
prediction. The results are constructed with M p = 40 and the weight of the regularization 
function 7 = 1. 



